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ABSTRACT 

The theoretical basis for well-posed marching of a Parabolic Navier-Stokes (PNS) 
computational technique for supersonic flow is discussed and examples given to verify 
the analysis. It is demonstrated that stable computations can be made even with very 
small steps in the marching direction. The method is applied to cones at large 
angles of attack in high Reynolds -number, supersonic flow. Streamline trajectories 
generated from the numerical solutions demonstrate the development of vortex 
structures on the lee side of the cone. 


INTRODUCTION 

The computation of steady supersonic viscous flows with one set of equations, 
valid in both the boundary layer and the inviscid supersonic region, has long been 
attempted to avoid the difficulties with matching conditions for separate inviscid 
and boundary-layer equations. Use of the unsteady form of the NS equations has been 
successful but inefficient, in terms of computer time and memory requirements. 
Therefore, many authors [1-8] have used steady marching or iterative methods to solve 
numerically the time- invariant NS equations. Single-pass methods have been used 
for supersonic flows, where the inviscid region is a well-posed initial-value 
problem, and where the viscous region is a boundary layer, which is itself amenable 
to a marching solution. The difficulty with solving the complete inviscid/viscous 
domain at once lies in the thin subsonic layer near the wall, where the pressure 
is determined by the solution. The absence of a downstream boundary condition makes 
the problem ill— posed. This combined problem is basically an Interaction of a 
boundary layer with an inviscid, supersonic stream. Lighthill [9] analyzed this 
problem and found solutions for the boundary-layer displacement thickness of the 
form 6 = a exp(kx). One can interpret this exponential growth as the onset of 
streamwise expansion or compression, which may lead to separation; these solutions 
have become known as '’departure** solutions corresponding to some undetermined 
downstream boundary condition. This paper describes an approach for avoiding the 
exponentially growing solution, and presents results for turbulent flow over a cone 
at large angles of attack. 


THE PNS EQUATIONS 

The governing equations are written in general curvilinear conservation-law 
form as follows [8]: 


^ {J- [5^E^ + ± {j"^n^(E - E^) + n„(F - F„) + n.(G - GJ]} 
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Here J is the Jacobian of the coordinate transormation; E, F, and G are the 
inviscid flux vectors; and and Gy are the viscous flux vectors — all written in 
terms of Cartesian velocity components. Note that the viscous terms are omitted 
from the ? derivative term; this is the fundamental parabolizing approximation. 
The ( )* terms are central to the discussion of parabolic marching, and are given 
by the column vectors 


* r 2 T 

E = Ipu, pu + cop, puv, puw, (pe^ + p)u} 


* 2 T 

F = {pv, puv, pv + cop, pvw, (pe^ + p)v} 


( 2 ) 


* 2 T 

G = {pw, puw, pvw, pw + cop , (pe^ + p)w} 


where 


pe^ = p / (y - 1 ) + 0 . 5 p(u^ + v^ + w^) 


Here u, v, and w are Cartesian velocity components, p is the pressure, p the 
density, and y ratio of specific heats. The P terms on the right side are 

given by 


= {0, (1 - co)p, 0, b, 0}^ , 

P2 = {0, 0, (1 - (0)p, 0, 0}^ , . 

P3 = {0, 0, 0, (1 - (0)p, 0}^ . 


( 3 ) 


Following Vigneron et al. [ 4 ], the parameter co is included with the pressure 
term in the momentum equations. It was shown that the equations can be made formally 
parabolic by utilizing an appropriate functional form for co. Considering the two- 
dimensional viscous subset of equations (1) , the mathematical character of the 
equations is governed by the eigenvalues of the differential system, leading to the 
condition 


U) < yM ^/[l + (y - 1 )M 

X X 


(A) 


If this condition is satisfied, then the equations remain parabolic even in 
the viscous region near the wall where the Mach number is less than unity; a similar 
form also governs the inviscid region, but is not pertinent to the present work. 

The purpose of this paper is to further demonstrate the effectiveness of this 
approach for interacting supersonic flows where the problem can be solved in a single 
sweep of the PNS equations. 

In previous work with the PNS method, a ^sublayer" approach was used, in which 
the pressure gradient term was completely removed in the subsonic region of flow. 

This is equivalent to a step function for cj, and this places greater significance 
on the P terms on the right side (RHS) of the equation. Usually, backward- 
difference approximations have been used for the pressure gradient term on the RHS. 
However, it can be shown that a backward difference leads to numerical instability 
for sufficiently small marching steps (which may be a manifestation of the departure 
behavior). When the o) function is used, the RHS can be set to zero, to a good 
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approximation, making the numerical computation stable, even for small steps. Thus, 
the (X) function and a zero RHS together permit the investigation of departure 
behavior, which requires very small streamwise steps. 


NUMERICAL METHODS 

The present computational results are obtained with the Beam-Warming factored 
algorithm in delta form. Details can be found in reference [8]. Full second— order 
accuracy can be maintained by use of central differences in n and and with a 
three-point backward difference in However, for the present results, the Euler 

implicit method is used for the marching direction, making the method first-order 
accurate in the ? direction. 


TRIPLE-DECK ANALYSIS 

For subsonic flows, that is, elliptic equations, a downstream boundary condition 
must be specified. When attempting to march such a system without specification 
of the downstream condition, there are an infinity of possible solutions that corre- 
spond to various attached or separated flows. When the viscous layer is thin, the 
departure solution is governed by the interaction of the viscous and inviscid regions 
near the wall [9]. Stewartson [10] further studied such flows using a triple-deck 
analysis consisting of a Llghthill sublayer, an inviscid shear region, and the 
inviscid external flow. This analysis establishes the order of magnitude of the 
interaction region. 

The pressure is constant across both the lower and middle decks, and the velocity 
in the middle deck corresponds with the attached velocity profile upstream of the 
interaction, but displaced by 6. Utilizing linear supersonic theory for the pressure 
at the top edge of the middle deck, one obtains the following momentum equation for 
the lower deck: 


3u Du _ Du 2, . 2 d^6 

“ 17 " " ■ ''u"u » 1 — 


dy 


dx 


(5) 


where the subscript w refers to wall, and e refers to the top edge of the middle 
deck. Introducing the stream function in (5), and utilizing the Lighthill form for 
the displacement thickness 6, yields an Airy equation. It is found from the 
solution that the triple-deck scale is given by 

Ak = 0.8272(M^ Ax = k Ax , (6) 

where x Is the nondimensional physical distance, Re is the Reynolds number, Cf 
is the skin friction coefficient, T is the temperature, and M is the Mach number. 
The pressure gradient that results is 

dp/dx = C exp(kx) . (7) 

The importance of this analysis to computations with the PNS equations is that 
it provides a test of the effectiveness of methods for damping the unwanted departure 
solutions. To Illustrate, we consider a flat-plate boundary layer. The solution 
shown in Figure 1 at x/L = 0.5 is obtained with the u parameter turned on, and 
then uj is set to 1.0, so that the full PNS equations are marched for a short 
distance. The solution immediately starts to "depart" until at x = 0.5055 the u 
factor is again turned on, and the solution decays toward the well-known boundary- 
layer solution. If the triple— deck scale had not been resolved by the step size 


3 



Ax, the departure would have appeared as an oscillatory divergence. We note here 
that other investigators report their methods are unable to suppress the departure 
behavior for small marching steps, while the present approach is quite effective. 

Figure 2 demonstrates that the computed departure solutions (oj = 1) agree with 
Lighthill’s linear analysis. This figure shows the magnitude of the computed 
pressure gradient as a function of the triple-deck distance scale x = kx. The 
initial slope of the computed departure agrees with the linear analysis until 
nonlinear effects take over. 


COMPUTATION OF TURBULENT FLOWS 

Turbulent computations are performed with the Reynolds averaged form of the 
PNS equations and an algebraic, or zero-equation, eddy viscosity model. Details are 
given in reference [5] and will not be repeated here. We note, however, that no 
special modifications were needed for the present calculations, even with a large 
crossflow separation on the leeward side. 

Figures 3 and 4 compare the circumferential variation of pressure and surface 
stream angle from the present computations with the experimental results of Rainbird 
[11]. Generally,, good agreement is obtained with experiment. The main differences 
are caused by an incorrect prediction of the location of separation, believed to be 
due to the inadequacy of the simple turbulence model used. These results were 
obtained on a CDC 7600 computer with 50 grid points between the body and the shock, 
and with 47 unequally spaced meridian planes. 


FLOW-FIELD SIMULATION 

The simulation of turbulent vortex structures is achieved by tracing particle 
paths using the velocity field obtained from the PNS solution. Initial positions 
for the tracer particles are specified at field locations near the cone surface, 
and x/L = 0.1. The particle paths are then determined from a simple Euler predictor 
corrector, finite-difference scheme. The velocity field is taken at x = L and is 
assumed invariant with x for the purposes of this simulation. 

The flow simulation is shown in Figures 5a and 5b. Three distinct vortices are 
observed — one emanating from the primary separation, the second from the secondary 
separation of the flow coming down the leeward plane of symmetry, and a third below 
the primary vortex and having the same rotation as the primary vortex. This flow 
structure is suggestive of a global type of mixing in large turbulent structures, 
which is significant to the development of turbulence models. 
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FLAT PLATE Moo = 3 Re = 50X106 



Fig. 1. Influence of "Omega" parameter on 
departure and recovery. 



Fig. 2. Departure solution for a flat plate; 
= 3, Re = 50x10®, T„ = 100 K. 
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